
# A1_Create_Baseline
#==============================================================================

# Description: This code collects the demand estimates and:
# - Extracts the parameter vector that minimized the objective 
# - If the code for a product category did not converge for the full sample, the
#   result for the specification without small shares is used

rm(list=ls())
setwd("D:/data_replication")   

starts <- 5
KK = 1060


# Import Demand Estimation Results for the full sample and Export to csv 
#===============================================================================

for (k in 1:KK) { 
objective_MPEC_starts = matrix(0,1,starts)
status_MPEC_starts = matrix(99,1,starts)             
filename <- sprintf("robustness/PC6/2_demand_estimation/2_results_full/resultsr_%s_%s.RData",k ,1) 
if (file.exists(filename)) {
load(filename)
X_MPEC_starts = X_MPEC 
objective_MPEC_starts[1,1] = objective_MPEC
status_MPEC_starts[1,1] = status_MPEC 
for (reps in 2:starts) {  
  filename <- sprintf("robustness/PC6/2_demand_estimation/2_results_full/resultsr_%s_%s.RData",k ,reps )
  X_MPEC <- matrix(0, nrow(X_MPEC_starts), 1)
  if (file.exists(filename)) {
  load(filename)    
  objective_MPEC_starts[1,reps] = objective_MPEC
  status_MPEC_starts[1,reps] = status_MPEC 
  }
  X_MPEC_starts <- cbind(X_MPEC_starts, X_MPEC)
}
X_MPEC = X_MPEC_starts
objective_MPEC = objective_MPEC_starts
status_MPEC = status_MPEC_starts
filename_X_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/X_MPEC_%s.csv", k)
write.csv(X_MPEC, file = filename_X_MPEC)
filename_objective_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/objective_MPEC_%s.csv", k)
write.csv(objective_MPEC, file = filename_objective_MPEC)
filename_status_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/status_MPEC_%s.csv", k)
write.csv(status_MPEC, file = filename_status_MPEC)
filename_data = sprintf("robustness/PC6/2_demand_estimation/3_baseline/data_%s.csv", k)
write.csv(data, file = filename_data)
filename_FE = sprintf("robustness/PC6/2_demand_estimation/3_baseline/FE_%s.csv", k)
write.csv(FE, file = filename_FE)
filename_v_save = sprintf("robustness/PC6/2_demand_estimation/3_baseline/v_save_%s.csv", k)
write.csv(v_save, file = filename_v_save)
filename_Y_save = sprintf("robustness/PC6/2_demand_estimation/3_baseline/Y_save_%s.csv", k)
write.csv(Y_save, file = filename_Y_save)
}
}

write.csv(K, file = "robustness/PC6/2_demand_estimation/3_baseline/K.csv")                                      # Save K to File
write.csv(K_A, file = "robustness/PC6/2_demand_estimation/3_baseline/K_A.csv")                                  # Save K_A to File


load("robustness/PC6/2_demand_estimation/codes_server/income_mu_sigma.RData")                 # Import Demographics
names(income_mu_sigma) <- c("declarant", "year", "mu", "sigma")
write.csv(income_mu_sigma, file = "robustness/PC6/2_demand_estimation/3_baseline/income_mu_sigma.csv", row.names=F)        # Export Demographics to csv

#===============================================================================


# Determine Products that converged 
#===============================================================================

converged <- matrix(0,KK,1)  
for (k in 1:KK) { 
  status_MPEC_starts = matrix(99,1,starts)  
  filename <- sprintf("robustness/PC6/2_demand_estimation/2_results_full/resultsr_%s_%s.RData",k ,1)  
  if (file.exists(filename)) {
  load(filename)
  status_MPEC_starts[1,1] = status_MPEC 
  for (reps in 2:starts) {  
    filename <- sprintf("robustness/PC6/2_demand_estimation/2_results_full/resultsr_%s_%s.RData",k ,reps )
    if (file.exists(filename)) {
    load(filename)
    status_MPEC_starts[1,reps] = status_MPEC 
    }
  }
  status_MPEC = status_MPEC_starts
  converged[k,1] = any(status_MPEC == 0)                                                                  # Test if any iteration converged (status == 0)
  }
}

#===============================================================================


# Replace results with those for the smaller sample if converged == 0 
#===============================================================================

converged_final <- converged

for (k in 1:KK) { 
  if (converged[k,1] == 0) {
  objective_MPEC_starts = matrix(0,1,starts)
  status_MPEC_starts = matrix(0,1,starts)  
  filename <- sprintf("robustness/PC6/2_demand_estimation/1_results_subsample/resultsr_%s_%s.RData",k ,1)  
  if (file.exists(filename)) {
  load(filename)
  X_MPEC_starts = X_MPEC 
  objective_MPEC_starts[1,1] = objective_MPEC
  status_MPEC_starts[1,1] = status_MPEC 
  for (reps in 2:starts) {  
    filename <- sprintf("robustness/PC6/2_demand_estimation/1_results_subsample/resultsr_%s_%s.RData",k ,reps) 
    if (file.exists(filename)) {
    load(filename)
    X_MPEC_starts <- cbind(X_MPEC_starts, X_MPEC)
    objective_MPEC_starts[1,reps] = objective_MPEC
    status_MPEC_starts[1,reps] = status_MPEC 
    }
  }
  X_MPEC = X_MPEC_starts
  objective_MPEC = objective_MPEC_starts
  status_MPEC = status_MPEC_starts
  filename_X_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/X_MPEC_%s.csv", k)
  write.csv(X_MPEC, file = filename_X_MPEC)
  filename_objective_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/objective_MPEC_%s.csv", k)
  write.csv(objective_MPEC, file = filename_objective_MPEC)
  filename_status_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/status_MPEC_%s.csv", k)
  write.csv(status_MPEC, file = filename_status_MPEC)
  filename_data = sprintf("robustness/PC6/2_demand_estimation/3_baseline/data_%s.csv", k)
  write.csv(data, file = filename_data)
  filename_FE = sprintf("robustness/PC6/2_demand_estimation/3_baseline/FE_%s.csv", k)
  write.csv(FE, file = filename_FE)
  filename_v_save = sprintf("robustness/PC6/2_demand_estimation/3_baseline/v_save_%s.csv", k)
  write.csv(v_save, file = filename_v_save)
  filename_Y_save = sprintf("robustness/PC6/2_demand_estimation/3_baseline/Y_save_%s.csv", k)
  write.csv(Y_save, file = filename_Y_save)
  converged_final[k,1] <- any(status_MPEC == 0) 
  }
  }
}

#===============================================================================


# Renumber product categories to account for those without convergence
#===============================================================================

pc8plus_i_k <- matrix(0, KK, 2)

i_k <- 0
for (k in 1:KK) { 
  if (converged_final[k,1] == 1) {
    i_k = i_k + 1
    pc8plus_i_k[k,1] = k
    pc8plus_i_k[k,2] = i_k
    filename_X_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/X_MPEC_%s.csv", k)
    filename_X_MPEC2 = sprintf("robustness/PC6/2_demand_estimation/3_baseline/X_MPEC_%s.csv", i_k)
    file.rename(filename_X_MPEC, filename_X_MPEC2)
    filename_objective_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/objective_MPEC_%s.csv", k)
    filename_objective_MPEC2 = sprintf("robustness/PC6/2_demand_estimation/3_baseline/objective_MPEC_%s.csv", i_k)
    file.rename(filename_objective_MPEC, filename_objective_MPEC2)
    filename_status_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/status_MPEC_%s.csv", k)
    filename_status_MPEC2 = sprintf("robustness/PC6/2_demand_estimation/3_baseline/status_MPEC_%s.csv", i_k)
    file.rename(filename_status_MPEC, filename_status_MPEC2)
    filename_data = sprintf("robustness/PC6/2_demand_estimation/3_baseline/data_%s.csv", k)
    filename_data2 = sprintf("robustness/PC6/2_demand_estimation/3_baseline/data_%s.csv", i_k)
    file.rename(filename_data, filename_data2)
    filename_FE = sprintf("robustness/PC6/2_demand_estimation/3_baseline/FE_%s.csv", k)
    filename_FE2 = sprintf("robustness/PC6/2_demand_estimation/3_baseline/FE_%s.csv", i_k)
    file.rename(filename_FE, filename_FE2)
    filename_v_save = sprintf("robustness/PC6/2_demand_estimation/3_baseline/v_save_%s.csv", k)
    filename_v_save2 = sprintf("robustness/PC6/2_demand_estimation/3_baseline/v_save_%s.csv", i_k)
    file.rename(filename_v_save, filename_v_save2)
    filename_Y_save = sprintf("robustness/PC6/2_demand_estimation/3_baseline/Y_save_%s.csv", k)
    filename_Y_save2 = sprintf("robustness/PC6/2_demand_estimation/3_baseline/Y_save_%s.csv", i_k)
    file.rename(filename_Y_save, filename_Y_save2)
  }
}
 
write.csv(file = 'robustness/PC6/2_demand_estimation/pc8plus_i_k.csv', pc8plus_i_k)


non_converged <- c(1016,1017,1037)

for (i in non_converged) {
  filename_X_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/X_MPEC_%s.csv", i)
  file.remove(filename_X_MPEC) 
  filename_objective_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/objective_MPEC_%s.csv", i)
  file.remove(filename_objective_MPEC)
  filename_status_MPEC = sprintf("robustness/PC6/2_demand_estimation/3_baseline/status_MPEC_%s.csv", i)
  file.remove(filename_status_MPEC)
  filename_data = sprintf("robustness/PC6/2_demand_estimation/3_baseline/data_%s.csv", i)
  file.remove(filename_data)
  filename_FE = sprintf("robustness/PC6/2_demand_estimation/3_baseline/FE_%s.csv", i)
  file.remove(filename_FE)
  filename_v_save = sprintf("robustness/PC6/2_demand_estimation/3_baseline/v_save_%s.csv", i)
  file.remove(filename_v_save)
  filename_Y_save = sprintf("robustness/PC6/2_demand_estimation/3_baseline/Y_save_%s.csv", i)
  file.remove(filename_Y_save)
}

